Paso 3.6

Análisis de clases latentes: modelos seleccionados, sin pueblo originario y sin año y recategorización de edad mujer y semana gestacional hito 1, caracterización de clases y medidas de ajuste (glca)

Show code

script src = "https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js"

Show code

 $(document).ready(function() {
    $('body').prepend('<div class=\"zoomDiv\"><img src=\"\" class=\"zoomImg\"></div>');
    // onClick function for all plots (img's)
    $('img:not(.zoomImg)').click(function() {
      $('.zoomImg').attr('src', $(this).attr('src')).css({width: '100%'});
      $('.zoomDiv').css({opacity: '1', width: 'auto', border: '1px solid white', borderRadius: '5px', position: 'fixed', top: '50%', left: '50%', marginRight: '-50%', transform: 'translate(-50%, -50%)', boxShadow: '0px 0px 50px #888888', zIndex: '50', overflow: 'auto', maxHeight: '100%'});
    });
    // onClick function for zoomImg
    $('img.zoomImg').click(function() {
      $('.zoomDiv').css({opacity: '0', width: '0%'}); 
    });
  });
  

Show code

<script src="hideOutput.js"></script> 

Show code

$(document).ready(function() {    
    $chunks = $('.fold');    
    $chunks.each(function () {      // add button to source code chunks     
    if ( $(this).hasClass('s') ) {       
        $('pre.r', this).prepend("<div class=\"showopt\">Show Source</div><br style=\"line-height:22px;\"/>");
            $('pre.r', this).children('code').attr('class', 'folded');     
            }      // add button to output chunks     
        if ( $(this).hasClass('o') ) {       
            $('pre:not(.r)', this).has('code').prepend("<div class=\"showopt\">Show Output</div><br style=\"line-height:22px;\"/>");       
            $('pre:not(.r)', this).children('code:not(r)').addClass('folded');        // add button to plots       
            $(this).find('img').wrap('<pre class=\"plot\"></pre>');       
            $('pre.plot', this).prepend("<div class=\"showopt\">Show Plot</div><br style=\"line-height:22px;\"/>");       
            $('pre.plot', this).children('img').addClass('folded');      
            }   
});    // hide all chunks when document is loaded   
    $('.folded').css('display', 'none')    // function to toggle the visibility   
    $('.showopt').click(function() {     
            var label = $(this).html();     
            if (label.indexOf("Show") >= 0) {       
                $(this).html(label.replace("Show", "Hide"));     
            } else {
              $(this).html(label.replace("Hide", "Show"));     
            }     
    $(this).siblings('code, img').slideToggle('fast', 'swing');   
    }); 
}); 

Cargamos los datos

Show code

rm(list = ls());gc()

          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  653032 34.9    1332051 71.2   901241 48.2
Vcells 1172494  9.0    8388608 64.0  2075142 15.9

Show code

#cargar glca con y sin resultado distal
load("data2_lca2_adj216_alt_sin_po_ano_2023_05_14.RData")

Cargamos los paquetes

Show code

knitr::opts_chunk$set(echo = TRUE)

if(!require(poLCA)){install.packages("poLCA")}
if(!require(poLCAParallel)){devtools::install_github("QMUL/poLCAParallel@package")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(parallel)){install.packages("parallel")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(tidyverse)){install.packages("tidyverse")}
try(if(!require(sjPlot)){install.packages("sjPlot")})
if(!require(emmeans)){install.packages("emmeans")}
if(!require(nnet)){install.packages("nnet")}
if(!require(here)){install.packages("here")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(progress)){install.packages("progress")}
if(!require(caret)){install.packages("caret")}
if(!require(rpart)){install.packages("rpart")}
if(!require(rpart.plot)){install.packages("rpart.plot")}
if(!require(partykit)){install.packages("partykit")}
if(!require(randomForest)){install.packages("randomForest")}
if(!require(ggcorrplot)){install.packages("ggcorrplot")}
if(!require(polycor)){install.packages("polycor")}
if(!require(tableone)){install.packages("tableone")}
if(!require(broom)){install.packages("broom")}
if(!require(plotly)){install.packages("plotly")}
if(!require(rsvg)){install.packages("rsvg")}
if(!require(DiagrammeRsvg)){install.packages("DiagrammeRsvg")}
if(!require(effects)){install.packages("effects")}
if(!require(glca)){install.packages("glca")}

Seleccionar modelos finales

Sin resultado distal

Show code

best_model_glca<- eval(parse(text = paste0("lca",best_model))) 
summary(best_model_glca)

Call:
glca(formula = f_preds, data = mydata_preds3, nclass = 5, n.init = 500, 
    decreasing = T, testiter = 500, maxiter = 10000, seed = seed, 
    verbose = FALSE)

Manifest items : CAUSAL EDAD_MUJER_REC PAIS_ORIGEN_REC HITO1_EDAD_GEST_SEM_REC MACROZONA PREV_TRAMO_REC 

Categories for manifest items :
                        Y = 1 Y = 2 Y = 3 Y = 4 Y = 5 Y = 6
CAUSAL                      2     3     4                  
EDAD_MUJER_REC              1     2     3     4     5      
PAIS_ORIGEN_REC             1     2     3                  
HITO1_EDAD_GEST_SEM_REC     1     2     3     4            
MACROZONA                   1     2     3     4     5     6
PREV_TRAMO_REC              1     2     3     4     5      

Model : Latent class analysis 

Number of latent classes : 5 
Number of observations : 3789 
Number of parameters : 104 

log-likelihood : -22209.52 
     G-squared : 1632.457 
           AIC : 44627.03 
           BIC : 45275.98 

Marginal prevalences for latent classes :
Class 1 Class 2 Class 3 Class 4 Class 5 
0.13706 0.54102 0.12651 0.04249 0.15292 

Class prevalences by group :
    Class 1 Class 2 Class 3 Class 4 Class 5
ALL 0.13706 0.54102 0.12651 0.04249 0.15292

Item-response probabilities :
CAUSAL 
         Y = 1  Y = 2  Y = 3
Class 1 0.4647 0.5353 0.0000
Class 2 0.4038 0.5962 0.0000
Class 3 0.1928 0.8072 0.0000
Class 4 0.0233 0.0000 0.9767
Class 5 0.0098 0.0000 0.9902
EDAD_MUJER_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5
Class 1 0.0043 0.0085 0.1703 0.5025 0.3144
Class 2 0.0072 0.0180 0.2106 0.4571 0.3071
Class 3 0.0000 0.0000 0.0405 0.4255 0.5340
Class 4 0.0000 0.1649 0.2335 0.4515 0.1501
Class 5 0.0017 0.3472 0.2463 0.2958 0.1089
PAIS_ORIGEN_REC 
         Y = 1  Y = 2  Y = 3
Class 1 0.0339 0.0431 0.9230
Class 2 0.0002 0.9998 0.0000
Class 3 0.0000 0.9252 0.0748
Class 4 0.0000 0.0000 1.0000
Class 5 0.0000 0.9935 0.0065
HITO1_EDAD_GEST_SEM_REC 
         Y = 1  Y = 2  Y = 3  Y = 4
Class 1 0.0230 0.1163 0.7339 0.1268
Class 2 0.0291 0.1767 0.6594 0.1348
Class 3 0.0150 0.3653 0.5698 0.0499
Class 4 0.0191 0.9809 0.0000 0.0000
Class 5 0.0089 0.9877 0.0035 0.0000
MACROZONA 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0059 0.5683 0.1051 0.0594 0.2083 0.0530
Class 2 0.0022 0.3208 0.1932 0.2004 0.0935 0.1899
Class 3 0.0000 0.7147 0.0922 0.0556 0.0698 0.0677
Class 4 0.0062 0.5563 0.0870 0.0166 0.2960 0.0379
Class 5 0.0042 0.3849 0.1747 0.1452 0.0862 0.2048
PREV_TRAMO_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5
Class 1 0.0123 0.0139 0.6086 0.2871 0.0780
Class 2 0.0012 0.0251 0.6227 0.3474 0.0036
Class 3 0.0000 0.8096 0.0000 0.1839 0.0065
Class 4 0.0253 0.0062 0.5225 0.1669 0.2791
Class 5 0.0017 0.0695 0.7236 0.1998 0.0053

Show code

#https://rdrr.io/cran/glca/src/R/summary.glca.R
#Class prevalence by group mean(best_model_glca_w_distal_outcome$posterior$ALL$`Class 1`)
#Item-response probabilities (most likely response)  print(sapply(best_model_glca_w_distal_outcome$param$rho[[1]],
#print(sapply(best_model_glca_w_distal_outcome$param$rho$ALL, function(m) apply(m, 1, which.max)))
#function(m) apply(m, 1, which.max)))

Show code

#https://rdrr.io/cran/glca/src/R/plot.glca.R
plot(eval(parse(text = paste0("lca",best_model))), ask=F)

Figure 1: Selected Model

Figure 2: Selected Model

Figure 3: Selected Model

Figure 4: Selected Model

Figure 5: Selected Model

Figure 6: Selected Model

Figure 7: Selected Model

Show code

rho_glca<- 
do.call("bind_rows",best_model_glca$param$rho$ALL) %>% 
  t() %>% 
  round(2) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  magrittr::set_colnames(c("variables", paste0("Class",1:length(best_model_glca$param$gamma)))) %>% 
  tidyr::separate(variables, into=c("var", "prob"), sep=".Y =")

lcmodel_glca <- reshape2::melt(rho_glca, level=2) %>% dplyr::rename("class"="variable")

traductor_cats <- readxl::read_excel("tabla12_corr.xlsx") %>% 
  dplyr::mutate(level=readr::parse_double(level)) %>% 
  dplyr::mutate(level= dplyr::case_when(grepl("CAUSAL",Name)~ level-1,T~level)) %>% 
  dplyr::mutate(level= dplyr::case_when(grepl("AÑO",Name)~ level-1,T~level)) %>% 
  dplyr::mutate(charactersitic=gsub(" \\(%\\)", "", Name))


lcmodel_glca<- lcmodel_glca %>% 
  dplyr::mutate(pr=as.numeric(gsub("[^0-9.]+", "", prob))) %>% 
  dplyr::left_join(traductor_cats[,c("charactersitic", "level", "CATEGORIA")], by= c("var"="charactersitic", "pr"="level"))  
  #dplyr::mutate(CATEGORIA= dplyr::case_when(var=="AÑO" & prob==" 1"~"Perdidos", T~CATEGORIA))

lcmodel_glca$text_label<-paste0("Categoria:",lcmodel_glca$CATEGORIA,"<br>%: ",scales::percent(lcmodel_glca$value))

zp3 <- ggplot(lcmodel_glca,aes(x = var, y = value, fill = factor(pr), label=text_label))
zp3 <- zp3 + geom_bar(stat = "identity", position = "stack")
zp3 <- zp3 + facet_grid(class ~ .) 
zp3 <- zp3 + scale_fill_brewer(type="seq", palette="Greys", na.value = "white") +theme_bw()
zp3 <- zp3 + labs(y = "Porcentaje de probabilidad de respuesta", 
                  x = "",
                  fill ="Cateorías de\nRespuesta")
zp3 <- zp3 + theme( axis.text.y=element_blank(),
                    axis.ticks.y=element_blank(),                    
                    panel.grid.major.y=element_blank())
zp3 <- zp3 + guides(fill = guide_legend(reverse=TRUE))
zp3 <- zp3 + theme(axis.text.x = element_text(angle = 30, hjust = 1))
#print(zp1)

ggplotly(zp3, tooltip = c("text_label"))%>% layout(xaxis= list(showticklabels = T),height=600, width=800)

Figure 8: Selected Model

Show code

ggsave("_fig3_LCA_distribuciones_glca_sin_po_ano.png",zp3, dpi= 600)

lcmodel_glca %>% rio::export("variables_probabilities_in_category_glca_sin_po_ano.xlsx")

Con resultado distal

Show code

best_model_glca_w_distal_outcome<- eval(parse(text = paste0("lca2",best_model2)))

summary(best_model_glca_w_distal_outcome)

Call:
glca(formula = f_adj, data = mydata_preds3, nclass = 5, n.init = 500, 
    decreasing = T, testiter = 500, maxiter = 10000, seed = seed, 
    verbose = FALSE)

Manifest items : CAUSAL EDAD_MUJER_REC PAIS_ORIGEN_REC HITO1_EDAD_GEST_SEM_REC MACROZONA PREV_TRAMO_REC 
Covariates (Level 1) : outcome 

Categories for manifest items :
                        Y = 1 Y = 2 Y = 3 Y = 4 Y = 5 Y = 6
CAUSAL                      2     3     4                  
EDAD_MUJER_REC              1     2     3     4     5      
PAIS_ORIGEN_REC             1     2     3                  
HITO1_EDAD_GEST_SEM_REC     1     2     3     4            
MACROZONA                   1     2     3     4     5     6
PREV_TRAMO_REC              1     2     3     4     5      

Model : Latent class analysis 

Number of latent classes : 5 
Number of observations : 3789 
Number of parameters : 108 

log-likelihood : -22149.63 
     G-squared : 2421.673 
           AIC : 44515.27 
           BIC : 45189.17 

Marginal prevalences for latent classes :
Class 1 Class 2 Class 3 Class 4 Class 5 
0.12504 0.54687 0.13337 0.04254 0.15218 

Class prevalences by group :
    Class 1 Class 2 Class 3 Class 4 Class 5
ALL 0.12504 0.54687 0.13337 0.04254 0.15218

Logistic regression coefficients :
            Class 1/5 Class 2/5 Class 3/5 Class 4/5
(Intercept)    0.4116    2.2509   -0.6620   -1.5768
outcome1      -0.6863   -1.1317    0.5664    0.3254
Item-response probabilities :
CAUSAL 
         Y = 1  Y = 2  Y = 3
Class 1 0.4626 0.5374 0.0000
Class 2 0.4068 0.5932 0.0000
Class 3 0.2021 0.7979 0.0000
Class 4 0.0245 0.0000 0.9755
Class 5 0.0049 0.0000 0.9951
EDAD_MUJER_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5
Class 1 0.0000 0.0089 0.1746 0.5022 0.3142
Class 2 0.0082 0.0181 0.2114 0.4563 0.3060
Class 3 0.0000 0.0000 0.0392 0.4345 0.5262
Class 4 0.0000 0.1647 0.2332 0.4514 0.1507
Class 5 0.0017 0.3483 0.2460 0.2951 0.1089
PAIS_ORIGEN_REC 
         Y = 1  Y = 2  Y = 3
Class 1 0.0121 0.0001 0.9878
Class 2 0.0059 0.9941 0.0000
Class 3 0.0000 0.9069 0.0931
Class 4 0.0000 0.0007 0.9993
Class 5 0.0000 0.9932 0.0068
HITO1_EDAD_GEST_SEM_REC 
         Y = 1  Y = 2  Y = 3  Y = 4
Class 1 0.0204 0.1142 0.7360 0.1294
Class 2 0.0290 0.1733 0.6603 0.1374
Class 3 0.0181 0.3704 0.5720 0.0396
Class 4 0.0192 0.9808 0.0000 0.0000
Class 5 0.0088 0.9877 0.0035 0.0000
MACROZONA 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5  Y = 6
Class 1 0.0000 0.5660 0.1060 0.0623 0.2147 0.0510
Class 2 0.0038 0.3205 0.1918 0.1998 0.0949 0.1893
Class 3 0.0000 0.7205 0.0937 0.0498 0.0700 0.0660
Class 4 0.0062 0.5563 0.0870 0.0165 0.2956 0.0384
Class 5 0.0038 0.3850 0.1751 0.1452 0.0862 0.2047
PREV_TRAMO_REC 
         Y = 1  Y = 2  Y = 3  Y = 4  Y = 5
Class 1 0.0098 0.0007 0.6260 0.2819 0.0817
Class 2 0.0020 0.0305 0.6254 0.3384 0.0037
Class 3 0.0000 0.7586 0.0043 0.2281 0.0091
Class 4 0.0256 0.0061 0.5214 0.1672 0.2796
Class 5 0.0018 0.0696 0.7237 0.1997 0.0053

Show code

rho_glca_adj<- 
do.call("bind_rows",best_model_glca_w_distal_outcome$param$rho$ALL) %>% 
  t() %>% 
  round(2) %>% 
  data.table::data.table(keep.rownames = T) %>% 
  magrittr::set_colnames(c("variables", paste0("Class",1:dim(best_model_glca_w_distal_outcome$param$gamma[[1]])[[2]]))) %>% 
  tidyr::separate(variables, into=c("var", "prob"), sep=".Y =")

lcmodel_glca_adj <- reshape2::melt(rho_glca_adj, level=2) %>% dplyr::rename("class"="variable")

lcmodel_glca_adj<- lcmodel_glca_adj %>% 
  dplyr::mutate(pr=as.numeric(gsub("[^0-9.]+", "", prob))) %>% 
  dplyr::left_join(traductor_cats[,c("charactersitic", "level", "CATEGORIA")], by= c("var"="charactersitic", "pr"="level")) %>% 
  dplyr::mutate(CATEGORIA= dplyr::case_when(var=="AÑO" & prob==" 1"~"Perdidos", T~CATEGORIA))

lcmodel_glca_adj$text_label<-paste0("Categoria:",lcmodel_glca_adj$CATEGORIA,"<br>%: ",scales::percent(lcmodel_glca_adj$value))

zp4 <- ggplot(lcmodel_glca_adj,aes(x = var, y = value, fill = factor(pr), label=text_label))
zp4 <- zp4 + geom_bar(stat = "identity", position = "stack")
zp4 <- zp4 + facet_grid(class ~ .) 
zp4 <- zp4 + scale_fill_brewer(type="seq", palette="Greys", na.value = "white") +theme_bw()
zp4 <- zp4 + labs(y = "Porcentaje de probabilidad de respuesta", 
                  x = "",
                  fill ="Cateorías de\nRespuesta")
zp4 <- zp4 + theme( axis.text.y=element_blank(),
                    axis.ticks.y=element_blank(),                    
                    panel.grid.major.y=element_blank())
zp4 <- zp4 + guides(fill = guide_legend(reverse=TRUE))
zp4 <- zp4 + theme(axis.text.x = element_text(angle = 30, hjust = 1))
#print(zp1)

ggplotly(zp4, tooltip = c("text_label"))%>% layout(xaxis= list(showticklabels = T),height=600, width=800)

Figure 9: Selected Model

Show code

ggsave("_fig4_LCA_distribuciones_glca_adj_sin_po_ano.png",zp4, dpi= 600)

lcmodel_glca %>% rio::export("variables_probabilities_in_category_glca_adj_sin_po_ano.xlsx")

Show code

#https://rdrr.io/cran/glca/src/R/plot.glca.R
plot(eval(parse(text = paste0("lca2",best_model2))), ask=F)

Figure 10: Selected Model

Figure 11: Selected Model

Figure 12: Selected Model

Figure 13: Selected Model

Figure 14: Selected Model

Figure 15: Selected Model

Figure 16: Selected Model

Show code

glca_gam_dist_outcome<- best_model_glca_w_distal_outcome$param$gamma[[1]]
glca_gam_dist_outcome<-colnames(paste0("Class",1:length(best_model_glca_w_distal_outcome$param$gamma)))
#conditional probabilities
#Pr(B1=1|Class 3)
posteriors_glca_adj <- data.frame(best_model_glca_w_distal_outcome$posterior, 
                             predclass=best_model_glca_w_distal_outcome$param$gamma) 

classification_table_adj <- plyr::ddply(posteriors, "predclass", function(x) colSums(x[,1:length(LCA_best_model_adj_mod$P)]))
clasification_errors_adj<- 1 - sum(diag(as.matrix(classification_table_adj[,2:(length(LCA_best_model_adj_mod$P)+1)]))) / sum(classification_table_adj[,2:(length(LCA_best_model_adj_mod$P)+1)]) 

warning(paste("Error de clasificación: ", round(clasification_errors_adj,2)))


entropy_alt <- function(p) sum(-p * log(p))
error_prior_adj <- entropy_alt(LCA_best_model_adj_mod$P) # Class proportions
error_post_adj <- mean(apply(LCA_best_model_adj_mod$posterior, 1, entropy_alt),na.rm=T)
R2_entropy_alt_adj <- (error_prior_adj - error_post_adj) / error_prior_adj
warning(paste("Entropía: ", round(R2_entropy_alt_adj,2)))


#https://stackoverflow.com/questions/72783185/entropy-calculation-gives-nan-is-applying-na-omit-a-valid-tweak
entropy.safe <- function (p) {
  if (any(p > 1 | p < 0)) stop("probability must be between 0 and 1")
  log.p <- numeric(length(p))
  safe <- p != 0
  log.p[safe] <- log(p[safe])
  sum(-p * log.p)
}

error_prior2_adj <- entropy.safe(LCA_best_model_adj_mod$P) # Class proportions
error_post2_adj <- mean(apply(LCA_best_model_adj_mod$posterior, 1, entropy.safe),na.rm=T)
R2_entropy_safe_adj <- (error_prior2_adj - error_post2_adj) / error_prior2_adj
warning(paste("Entropía segura: ", round(R2_entropy_safe,2)))

#https://gist.github.com/daob/c2b6d83815ddd57cde3cebfdc2c267b3
warning(paste("Entropía (solución Oberski): ", round(entropy.R2(LCA_best_model_adj_mod),2)))

#\#minimum average posterior robability of cluster membership (\>0.7), interpretability (classes are clearly distinguishable), and parsimony (each class has a sufficient sample size for further analysis; n≥50).

Show code

coef(eval(parse(text = paste0("lca2",best_model2))))

Class 1 / 5 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)    
(Intercept)     2.9257      1.0735      0.3718    2.887  0.003907 ** 
outcome1        0.2857     -1.2527      0.3355   -3.734  0.000191 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 2 / 5 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)  
(Intercept)     0.4006     -0.9149      0.5006   -1.828    0.0677 .
outcome1        0.7858     -0.2410      0.4769   -0.505    0.6133  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 3 / 5 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)    
(Intercept)    18.4092      2.9128      0.3577    8.142  5.26e-16 ***
outcome1        0.1830     -1.6981      0.3147   -5.396  7.26e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Class 4 / 5 :
            Odds Ratio Coefficient  Std. Error  t value  Pr(>|t|)  
(Intercept)     1.9386      0.6620      0.3775    1.754    0.0796 .
outcome1        0.5676     -0.5664      0.3464   -1.635    0.1021  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Show code

save.image("data2_lca3_glca_sin_po_ano.RData")

Comparar modelos

Show code

#Posterior mean class
round(apply(best_model_glca$posterior$ALL, 2,mean)*100,2)

Class 1 Class 2 Class 3 Class 4 Class 5 
  13.71   54.10   12.65    4.25   15.29 

Show code

#Classifying by posterior probs.
posterior_glca_07<-
best_model_glca$posterior$ALL %>% 
    dplyr::mutate_all(~ifelse(.>.7,1,0)) %>% 
  dplyr::mutate(final_07=dplyr::case_when(`Class 1`==1~1,`Class 2`==1~2, `Class 3`==1~3,`Class 4`==1~4, `Class 5`==1~5))
posterior_glca_07 %>% 
    rowwise() %>%
  mutate(count_ones = sum(c_across(starts_with("Class")) == 1)) %>%
  ungroup() %>% 
  janitor::tabyl(final_07,count_ones)

 final_07   0    1
        1   0  499
        2   0 2035
        3   0  389
        4   0  161
        5   0  571
       NA 134    0

Show code

posterior_glca_05<-
best_model_glca$posterior$ALL %>% 
    dplyr::mutate_all(~ifelse(.>.5,1,0)) %>% 
  dplyr::mutate(final_05=dplyr::case_when(`Class 1`==1~1,`Class 2`==1~2, `Class 3`==1~3,`Class 4`==1~4, `Class 5`==1~5))
posterior_glca_05 %>% 
  rowwise() %>%
  mutate(count_ones = sum(c_across(starts_with("Class")) == 1)) %>%
  ungroup() %>% 
  janitor::tabyl(final_05, count_ones)

 final_05 0    1
        1 0  502
        2 0 2106
        3 0  447
        4 0  161
        5 0  571
       NA 2    0

Show code

#comparar clasificaciones

load("data2_lca3_sin_po_ano_2023_05_14.RData")

posterior_polca_05<-
LCA_best_model_mod$posterior %>% 
  data.frame() %>% 
  dplyr::mutate_all(~ifelse(.>.5,1,0)) %>% 
  dplyr::mutate(final_05=dplyr::case_when(`X1`==1~1,`X2`==1~2, `X3`==1~3,`X4`==1~4, `X5`==1~5))

table(posterior_polca_05$final_05,LCA_best_model_mod$predclass)

   
       1    2    3    4    5
  1  161    0    0    0    0
  2    0  447    0    0    0
  3    0    0  502    0    0
  4    0    0    0 2106    0
  5    0    0    0    0  571

Show code

warning(paste0("Sujetos fuera de la clasificación >0.5: ",round(1-sum(diag(table(posterior_polca_05$final_05,LCA_best_model_mod$predclass, exclude=NULL)))/length(LCA_best_model_adj_mod$predclass),3)))

posterior_polca_07<-
LCA_best_model_mod$posterior %>% 
  data.frame() %>% 
  dplyr::mutate_all(~ifelse(.>.7,1,0)) %>% 
  dplyr::mutate(final_07=dplyr::case_when(`X1`==1~1,`X2`==1~2, `X3`==1~3,`X4`==1~4, `X5`==1~5))

table(posterior_polca_07$final_07,LCA_best_model_mod$predclass, exclude=NULL)

      
          1    2    3    4    5
  1     161    0    0    0    0
  2       0  389    0    0    0
  3       0    0  499    0    0
  4       0    0    0 2035    0
  5       0    0    0    0  571
  <NA>    0   58    3   73    0

Show code

warning(paste0("Sujetos fuera de la clasificación >0.7: ",round(1-sum(diag(table(posterior_polca_07$final_07,LCA_best_model_mod$predclass, exclude=NULL)))/length(LCA_best_model_mod$predclass),2)))

LCA_best_model_mod$posterior %>% 
  data.frame() %>% 
  dplyr::mutate_all(~ifelse(.>.7,1,0)) %>% 
  dplyr::mutate(final_07=dplyr::case_when(`X1`==1~1,`X2`==1~2, `X3`==1~3,`X4`==1~4, `X5`==1~5))%>% 
  rowwise() %>%
  mutate(count_ones = sum(c_across(starts_with("Class")) == 1)) %>%
  ungroup() %>% 
  janitor::tabyl(final_07, count_ones)

 final_07    0
        1  161
        2  389
        3  499
        4 2035
        5  571
       NA  134

Show code

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
invisible("Cuán diferentes son unos modelos de otros??")
table(posterior_glca_07$final_07,posterior_polca_07$final_07, exclude=NULL)

      
          1    2    3    4    5 <NA>
  1       0    0  499    0    0    0
  2       0    0    0 2035    0    0
  3       0  389    0    0    0    0
  4     161    0    0    0    0    0
  5       0    0    0    0  571    0
  <NA>    0    0    0    0    0  134

Show code

table(posterior_glca_07$final_07,car::recode(posterior_polca_07$final_07,"3=1;4=2;2=3;1=4;NA=NA"),exclude=NULL)

      
          1    2    3    4    5 <NA>
  1     499    0    0    0    0    0
  2       0 2035    0    0    0    0
  3       0    0  389    0    0    0
  4       0    0    0  161    0    0
  5       0    0    0    0  571    0
  <NA>    0    0    0    0    0  134

Show code

sum(diag(table(posterior_glca_07$final_07,car::recode(posterior_polca_07$final_07,"3=1;4=2;2=3;1=4;NA=NA"))))/sum(table(posterior_glca_07$final_07,car::recode(posterior_polca_07$final_07,"3=1;4=2;2=3;1=4;NA=NA"),exclude=NULL))

[1] 0.9646345

Show code

sum(diag(table(posterior_glca_07$final_07,car::recode(posterior_polca_07$final_07,"3=1;4=2;2=3;1=4;NA=NA"))))/sum(table(posterior_glca_07$final_07,posterior_polca_07$final_07, exclude=NULL))

[1] 0.9646345

Show code

paste0("Perdidos en ambos: ", data.frame(table(posterior_glca_07$final_07,posterior_polca_07$final_07, exclude=NULL))[36,3])

[1] "Perdidos en ambos: 134"

Show code

table(posterior_glca_05$final_05,posterior_polca_05$final_05, exclude=NULL)

      
          1    2    3    4    5 <NA>
  1       0    0  502    0    0    0
  2       0    0    0 2106    0    0
  3       0  447    0    0    0    0
  4     161    0    0    0    0    0
  5       0    0    0    0  571    0
  <NA>    0    0    0    0    0    2

Show code

sum(diag(table(posterior_glca_05$final_05,car::recode(posterior_polca_05$final_05,"3=1;4=2;2=3;1=4;NA=NA"))))/sum(table(posterior_glca_05$final_05,posterior_polca_05$final_05, exclude=NULL))

[1] 0.9994722

Show code

paste0("Perdidos en ambos: ", data.frame(table(posterior_glca_05$final_05,posterior_polca_05$final_05, exclude=NULL))[36,3])

[1] "Perdidos en ambos: 2"

Show code

#best_model_glca_w_distal_outcome
#LCA_best_model_adj_mod

Show code

#Posterior mean class
round(apply(best_model_glca$posterior$ALL, 2,mean)*100,2)

Class 1 Class 2 Class 3 Class 4 Class 5 
  13.71   54.10   12.65    4.25   15.29 

Show code

#Classifying by posterior probs.
posterior_glca_07adj<-
best_model_glca_w_distal_outcome$posterior$ALL %>% 
    dplyr::mutate_all(~ifelse(.>.7,1,0)) %>% 
  dplyr::mutate(final_07=dplyr::case_when(`Class 1`==1~1,`Class 2`==1~2, `Class 3`==1~3,`Class 4`==1~4, `Class 5`==1~5))
posterior_glca_07adj %>% 
    rowwise() %>%
  mutate(count_ones = sum(c_across(starts_with("Class")) == 1)) %>%
  ungroup() %>% 
  janitor::tabyl(final_07,count_ones)

 final_07   0    1
        1   0  475
        2   0 2003
        3   0  382
        4   0  161
        5   0  571
       NA 197    0

Show code

posterior_glca_05adj<-
best_model_glca_w_distal_outcome$posterior$ALL %>% 
    dplyr::mutate_all(~ifelse(.>.5,1,0)) %>% 
  dplyr::mutate(final_05=dplyr::case_when(`Class 1`==1~1,`Class 2`==1~2, `Class 3`==1~3,`Class 4`==1~4, `Class 5`==1~5))
posterior_glca_05 %>% 
  rowwise() %>%
  mutate(count_ones = sum(c_across(starts_with("Class")) == 1)) %>%
  ungroup() %>% 
  janitor::tabyl(final_05, count_ones)

 final_05 0    1
        1 0  502
        2 0 2106
        3 0  447
        4 0  161
        5 0  571
       NA 2    0

Show code

#comparar clasificaciones

posterior_polca_05adj<-
LCA_best_model_adj_mod$posterior %>% 
  data.frame() %>% 
  dplyr::mutate_all(~ifelse(.>.5,1,0)) %>% 
  dplyr::mutate(final_05=dplyr::case_when(`X1`==1~1,`X2`==1~2, `X3`==1~3,`X4`==1~4, `X5`==1~5))

table(posterior_polca_05adj$final_05,LCA_best_model_adj_mod$predclass)

   
       1    2    3    4    5
  1  571    0    0    0    0
  2    0  488    0    0    0
  3    0    0 2116    0    0
  4    0    0    0  161    0
  5    0    0    0    0  453

Show code

warning(paste0("Sujetos fuera de la clasificación >0.5: ",round(1-sum(diag(table(posterior_polca_05adj$final_05,LCA_best_model_adj_mod$predclass, exclude=NULL)))/length(LCA_best_model_adj_mod$predclass),3)))

posterior_polca_07adj<-
LCA_best_model_adj_mod$posterior %>% 
  data.frame() %>% 
  dplyr::mutate_all(~ifelse(.>.7,1,0)) %>% 
  dplyr::mutate(final_07=dplyr::case_when(`X1`==1~1,`X2`==1~2, `X3`==1~3,`X4`==1~4, `X5`==1~5))

table(posterior_polca_07adj$final_07,LCA_best_model_adj_mod$predclass, exclude=NULL)

      
          1    2    3    4    5
  1     571    0    0    0    0
  2       0  475    0    0    0
  3       0    0 2003    0    0
  4       0    0    0  161    0
  5       0    0    0    0  382
  <NA>    0   13  113    0   71

Show code

warning(paste0("Sujetos fuera de la clasificación >0.7: ",round(1-sum(diag(table(posterior_polca_07adj$final_07,LCA_best_model_adj_mod$predclass)))/length(LCA_best_model_adj_mod$predclass),2)))

LCA_best_model_adj_mod$posterior %>% 
  data.frame() %>% 
  dplyr::mutate_all(~ifelse(.>.7,1,0)) %>% 
  dplyr::mutate(final_07=dplyr::case_when(`X1`==1~1,`X2`==1~2, `X3`==1~3,`X4`==1~4, `X5`==1~5))%>% 
  rowwise() %>%
  mutate(count_ones = sum(c_across(starts_with("Class")) == 1)) %>%
  ungroup() %>% 
  janitor::tabyl(final_07, count_ones)

 final_07    0
        1  571
        2  475
        3 2003
        4  161
        5  382
       NA  197

Show code

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
invisible("Cuán diferentes son unos modelos de otros??")
table(posterior_glca_07adj$final_07,posterior_polca_07adj$final_07, exclude=NULL)

      
          1    2    3    4    5 <NA>
  1       0  475    0    0    0    0
  2       0    0 2003    0    0    0
  3       0    0    0    0  382    0
  4       0    0    0  161    0    0
  5     571    0    0    0    0    0
  <NA>    0    0    0    0    0  197

Show code

table(posterior_glca_07adj$final_07,car::recode(posterior_polca_07adj$final_07,"1=5;2=1;3=2;4=4;5=3;NA=NA"),exclude=NULL)

      
          1    2    3    4    5 <NA>
  1     475    0    0    0    0    0
  2       0 2003    0    0    0    0
  3       0    0  382    0    0    0
  4       0    0    0  161    0    0
  5       0    0    0    0  571    0
  <NA>    0    0    0    0    0  197

Show code

sum(diag(table(posterior_glca_07adj$final_07,car::recode(posterior_polca_07adj$final_07,"1=5;2=1;3=2;4=4;5=3;NA=NA"))))/sum(table(posterior_glca_07adj$final_07,car::recode(posterior_polca_07adj$final_07,"1=5;2=1;3=2;4=4;5=3;NA=NA"),exclude=NULL))

[1] 0.9480074

Show code

paste0("Perdidos en ambos: ", data.frame(table(posterior_glca_07adj$final_07,posterior_polca_07adj$final_07, exclude=NULL))[36,3])

[1] "Perdidos en ambos: 197"

Show code

sum(diag(table(posterior_glca_07adj$final_07,car::recode(posterior_polca_07adj$final_07,"1=5;2=1;3=2;4=4;5=3;NA=NA"))))/sum(table(posterior_glca_07adj$final_07,posterior_polca_07adj$final_07, exclude=NULL))

[1] 0.9480074

Show code

table(posterior_glca_05adj$final_05,posterior_polca_05adj$final_05, exclude=NULL)

   
       1    2    3    4    5
  1    0  488    0    0    0
  2    0    0 2116    0    0
  3    0    0    0    0  453
  4    0    0    0  161    0
  5  571    0    0    0    0

Show code

sum(diag(table(posterior_glca_05adj$final_05,car::recode(posterior_polca_05adj$final_05,"1=5;2=1;3=2;4=4;5=3;NA=NA"))))/sum(table(posterior_glca_05adj$final_05,posterior_polca_05adj$final_05, exclude=NULL))

[1] 1

Show code

paste0("Perdidos en ambos: ", data.frame(table(posterior_glca_05adj$final_05,posterior_polca_05adj$final_05, exclude=NULL))[36,3])

[1] "Perdidos en ambos: NA"

Show code

invisible("ojo que no calzan porque son perdidos")

#best_model_glca_w_distal_outcome
#LCA_best_model_adj_mod

Show code

require(tidyverse)
sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
) %>% 
  DT::datatable(filter = 'top', colnames = c('Row number' =1,'Variable' = 2, 'Percentage'= 3),
              caption = htmltools::tags$caption(
        style = 'caption-side: top; text-align: left;',
        '', htmltools::em('Packages')),
      options=list(
initComplete = htmlwidgets::JS(
        "function(settings, json) {",
        "$(this.api().tables().body()).css({
            'font-family': 'Helvetica Neue',
            'font-size': '50%', 
            'code-inline-font-size': '15%', 
            'white-space': 'nowrap',
            'line-height': '0.75em',
            'min-height': '0.5em'
            });",#;
        "}")))